## Admissions-as-corrections reduce support for partisan misperceptions 
#and intended partisan media consumption (pilot anlyses)
##(Ms. No. 240015)--The Journal of Politics##
##Replication code for JOP Dataverse-pilot##
##Prepared by J.J. Fahey, 4/29/2025##

##The following libraries are necessary for analysis and data viz##
library(rio)
library(foreign)
library(tidyverse)
library(tidyr)
library(dplyr)
library(estimatr)
library(modelsummary)
library(stargazer)
library(lavaan)
library(pwr)
library(ggplot2)
library(nnet)
library(MNLpred)
library(scales)
library(knitr)

#Set working directory and import data

setwd("~/Dropbox/Research/Current Research/Admissions/JOP Reg Report/JOP Submission/JOP Data")

####Calculate differential attrition####
#Import data with incompletes

fox.pilot.fails <- import("Pilot_incompletes.xlsx")

#Sort into treatment and control conditions 
##Treatment Condition Assignment--admissions = 0 means control
#admissions = 1 is treatment##
fox.pilot.fails <- fox.pilot.fails %>%
  mutate(admissions = case_when (
    fox.pilot.fails$FL_22_DO_FL_23 == 1 ~ 0, 
    fox.pilot.fails$FL_22_DO_FL_24 == 1 ~ 1
  ))

##Dummy variable--whether respondent completed survey or not 
fox.pilot.fails <- fox.pilot.fails %>%
  mutate (complete = case_when (
    fox.pilot.fails$Progress == 100 ~ 1,
    fox.pilot.fails$Progress < 100 ~ 0
  ))

##Table of complete/incomplete by treatment (Appendix C.3)
table(fox.pilot.fails$admissions, fox.pilot.fails$complete)
t.test(complete ~ admissions, data = fox.pilot.fails)

####Calculating AV attention check failures####
##Import data with all observations

fox.pilot.avattn <- import("Pilot_all.xlsx")

fox.pilot.avattn <- fox.pilot.avattn %>%
  mutate(avattncheck = case_when(
    fox.pilot.avattn$avattn == 3 ~ 1,
    fox.pilot.avattn$avattn == 1 ~ 0,
    fox.pilot.avattn$avattn == 2 ~ 0,
    fox.pilot.avattn$avattn == 4 ~ 0))

##Calculation of number of respondents removed due to failing avattn check
#Appendix C.2 (column 1)
mean(fox.pilot.avattn$avattncheck, na.rm = TRUE)

##Note: prior to import, manually removed non-completes (Voluntary and those who failed screening questions) from dataset

####Import main (cleaned) data for analysis####
fox.pilot <- import("Pilot_cleaned.xlsx")

##Treatment Condition Assignment##
fox.pilot <- fox.pilot %>%
  mutate(admissions = case_when (
    fox.pilot$FL_22_DO == "FL_23" ~ 0, 
    fox.pilot$FL_22_DO == "FL_24" ~ 1
  ))

##Rudy vs. Sidney Assignment## 
fox.pilot <- fox.pilot %>%
  mutate(rudy = case_when (
    fox.pilot$FL_58_DO == "FoxVideoDobbs" ~ 1, 
    fox.pilot$FL_58_DO == "FoxVideoBartiromo" ~ 0,
    fox.pilot$FL_59_DO == "FoxVideoDobbs" ~ 1, 
    fox.pilot$FL_59_DO == "FoxVideoBartiromo" ~ 0,
  ))

#### Checks####

##Attn Check 1##
fox.pilot <- fox.pilot %>%
  mutate(attncheck1= case_when (
    fox.pilot$attn1 == 4 ~ 1, 
    fox.pilot$attn1 == 3 ~ 0, 
    fox.pilot$attn1 == 2 ~ 0, 
    fox.pilot$attn1 == 1 ~ 0, 
    fox.pilot$attn1 == 0 ~ 0
  ))

##Attn Check 2## 
fox.pilot <- fox.pilot %>%
  mutate(attncheck2= case_when (
    fox.pilot$attn2 == 1 ~ 1, 
    fox.pilot$attn2 == 0 ~ 0
  ))

##Attn Check 3##
fox.pilot <- fox.pilot %>%
  mutate(attncheck3= case_when (
    fox.pilot$attn3 == 1 ~ 1, 
    fox.pilot$attn3 == 0 ~ 0
  ))


#In order of discriminatory power, 3 is hardest, then 1, then 2#

##Dropping respondents who failed more than one attncheck 
fox.pilot <- fox.pilot %>%
  mutate(attntotal = attncheck1 + attncheck2 + attncheck3)

fox.pilot <- fox.pilot[fox.pilot$attntotal >= 2, ]

##Dropping respondents who indicate they are democrats 
fox.pilot <- fox.pilot[fox.pilot$ind1 == 0 | fox.pilot$ind1 == 2 | fox.pilot$pid == 1,]

####Pre-Analysis: Demographic Variable manipulation####

##Change age to numeric variable##

fox.pilot$age <- as.numeric(fox.pilot$age)
is.numeric(fox.pilot$age)

##Four point PID scale (because Dems are dropped)##
fox.pilot <- fox.pilot %>%
  mutate(pid.scale = case_when(
    fox.pilot$rep1 == 1~1, 
    fox.pilot$rep1 == 0~.666, 
    fox.pilot$ind1 == 0 ~ .333,
    fox.pilot$ind1 == 2 ~ 0))

##White dummy##

fox.pilot$white <- ifelse(grepl("^0$", fox.pilot$race), 1, 0)

##Hispanic/Latino dummy##
fox.pilot$hispanic <- ifelse(grepl("^0$", fox.pilot$hispanic), 1, 0)

##College dummmy## 

fox.pilot <- fox.pilot %>%
  mutate(college = case_when (
    fox.pilot$educ <= 4 ~ 0, 
    fox.pilot$educ >= 5 ~ 1
  ))

##Male Dummy## 
fox.pilot <- fox.pilot %>%
  mutate(male = case_when (
    fox.pilot$gender == 0 ~1, 
    fox.pilot$gender == 1 ~0, 
    fox.pilot$gender == 2 ~0))

##TV Station Dummies## 
##Fox##
fox.pilot <- fox.pilot %>%
  mutate(fox.viewer = ifelse(grepl("1", tvstations), 1, 0))

##OANN## 
fox.pilot <- fox.pilot %>%
  mutate(oann.viewer = ifelse(grepl("2", tvstations), 1, 0))

##CBS## 
fox.pilot <- fox.pilot %>%
  mutate(CBS.viewer = ifelse(grepl("3", tvstations), 1, 0))

##MSNBC## 
fox.pilot <- fox.pilot %>%
  mutate(MSNBC.viewer = ifelse(grepl("4", tvstations), 1, 0))

##CNN##
fox.pilot <- fox.pilot %>%
  mutate(CNN.viewer = ifelse(grepl("5", tvstations), 1, 0))

##ABC##
fox.pilot <- fox.pilot %>%
  mutate(ABC.viewer = ifelse(grepl("6", tvstations), 1, 0))

##Newsmax## 
fox.pilot <- fox.pilot %>%
  mutate(newsmax.viewer = ifelse(grepl("7", tvstations), 1, 0))

##Media trust## 
fox.pilot <- fox.pilot %>%
  mutate(media.trust = case_when (
    fox.pilot$mediatrust== 0 ~ 0,
    fox.pilot$mediatrust == 1 ~ .25,
    fox.pilot$mediatrust == 2 ~ .5,
    fox.pilot$mediatrust == 3 ~ .75,
    fox.pilot$mediatrust== 4 ~ 1))

##Calculating mean for summary statistics table (Appendix C.1) 

mean(fox.pilot$age, na.rm = TRUE)
mean(fox.pilot$college, na.rm = TRUE)
mean(fox.pilot$white, na.rm = TRUE)
mean(fox.pilot$male, na.rm = TRUE)

##Calculating Passage Rates for Attention Checks##

# Calculate means for each variable
means <- fox.pilot %>%
  summarise(
    Agree = mean(attncheck1, na.rm = TRUE),
    Color = mean(attncheck2, na.rm = TRUE),
    Matrix = mean(attncheck3, na.rm = TRUE)
  )

#Appendix C.2
kable(means, format = "latex", booktabs = TRUE, 
      caption = "Mean Scores for Attention Check Variables")


####Pre-Analysis: Mediator variables manipulation####
##Trust in Fox News##
fox.pilot <- fox.pilot %>%
  mutate(fox.trust = case_when(
    fox.pilot$matrixtrust_1 == 0 ~ 0,
    fox.pilot$matrixtrust_1 == 1 ~ .25,
    fox.pilot$matrixtrust_1 == 2 ~ .5,
    fox.pilot$matrixtrust_1 == 3 ~ .75,
    fox.pilot$matrixtrust_1 == 4 ~ 1))

##Partisan Identity Threat## 
fox.pilot <- fox.pilot %>%
  mutate(id.threat = case_when(
    fox.pilot$subjectivethreat == 0 ~ 0,
    fox.pilot$subjectivethreat == 1 ~ .25,
    fox.pilot$subjectivethreat == 2 ~ .5,
    fox.pilot$subjectivethreat == 3 ~ .75,
    fox.pilot$subjectivethreat == 4 ~ 1))

####Pre-Analysis: Dependent Variables Manipulation####
##Trump won##
fox.pilot <- fox.pilot %>%
  mutate(trump.won = case_when (
    fox.pilot$electionwinner == 1 ~ 1, 
    fox.pilot$electionwinner == 2 ~ 0, 
    fox.pilot$electionwinner == 0 ~ 0))

##Dominion changed votes## 
fox.pilot <- fox.pilot %>%
  mutate(dominion.changed = case_when (
    fox.pilot$voteschanged == 1 ~ 1, 
    fox.pilot$voteschanged == 2 ~ 0, 
    fox.pilot$voteschanged == 0 ~ 0))

##Immigrants Voted## 
fox.pilot <- fox.pilot %>%
  mutate(imm.voted = case_when (
    fox.pilot$immvoted == 1 ~ 1, 
    fox.pilot$immvoted == 2 ~ 0, 
    fox.pilot$immvoted == 0 ~ 0))

##venezuela## 
fox.pilot <- fox.pilot %>%
  mutate(venezuela.company= case_when (
    fox.pilot$venezuela == 0 ~ 1, 
    fox.pilot$venezuela == 2 ~ 0, 
    fox.pilot$venezuela == 1 ~ 0))

####Pre Analysis: Do Beliefs in Big Lie Form a Factor (CFA)####
##Appendix C.4

cfamodel2  <- ' h =~ NA*trump.won + dominion.changed + imm.voted + venezuela.company
h ~~ 1*h ' 

fit.cfa2 <- sem(cfamodel2, data = fox.pilot)
summary(fit.cfa2, fit.measures = TRUE, standardized = TRUE)

##Creating New Omnibus Variable## 

fox.pilot$big.lie <- (fox.pilot$trump.won + fox.pilot$dominion.changed + fox.pilot$imm.voted + 
  fox.pilot$venezuela.company)/4

####Main Analysis: Effect of Admissions as Corrections on Belief in Big Lie####

###Regression Model 1: Treatment on Support for Factor and Individual Variables (simple)###
##Appendix C.6
models_simple <- list(
  "Big Lie" = lm(big.lie ~ admissions, data = fox.pilot),
  "Trump Won" = lm(trump.won ~ admissions, data = fox.pilot),
  "Dominion Changed" = lm(dominion.changed ~ admissions, data = fox.pilot),
  "Immigrants Voted" = lm(imm.voted ~ admissions, data = fox.pilot), 
  "Venezuela" = lm(venezuela.company ~ admissions, data = fox.pilot))

modelsummary(models_simple,
             stars = TRUE,
             coef_rename = c("admissions" = "Admissions"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch, Nov. 25th-26th, 2023.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

###Regression Model 2: Treatment on Support for Factor and Individual Variables w/covariates###
##Appendix C.7
models_covariates <- list(
  "Big Lie"     = lm_lin(big.lie ~ admissions, ~ white + college + male + age + fox.viewer
                     + media.trust + pid.scale, data = fox.pilot),
  "Trump Won" = lm_lin(trump.won ~ admissions, ~ white + college + male + age + fox.viewer
                   + media.trust + pid.scale, data = fox.pilot),
  "Dominion Changed" = lm_lin(dominion.changed ~ admissions, ~ white + college + male + age + fox.viewer
                          + media.trust + pid.scale, data = fox.pilot),
  "Immigrants Voted" = lm_lin(imm.voted ~ admissions, ~  white + college + male + age + fox.viewer
                          + media.trust + pid.scale, data = fox.pilot), 
  "Venezuela Owns" = lm_lin(venezuela.company ~ admissions, ~ white + college + male + age + fox.viewer
                        + media.trust + pid.scale, data = fox.pilot))

modelsummary(models_covariates,
             stars = TRUE,
             coef_rename = c("admissions" = "Admissions", "white" = "white", 
                             "college" = "College educated", "male" = "male", 
                             "age" = "age", "fox.viewer" = "Fox viewer",
                             "media.trust" = "Media trust", "pid.scale" = "PID Scale"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch, Nov. 25th-26th, 2023.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

####Main Analysis: Effect of Treatment on Trust in Fox News and Subjective Identity Threat (hypothesized mediators)####

##Rename and manipulate trust in Fox News##
fox.pilot <- fox.pilot %>%
  mutate(fox.trust= case_when (
    fox.pilot$matrixtrust_1 == 5 ~ 1, 
    fox.pilot$matrixtrust_1 == 4 ~ .75, 
    fox.pilot$matrixtrust_1 == 3 ~ .5, 
    fox.pilot$matrixtrust_1 == 2 ~ .25, 
    fox.pilot$matrixtrust_1 == 1 ~ 0
  ))

##Rename and manipulate subjective identity threat## 
fox.pilot <- fox.pilot %>%
  mutate(id.threat= case_when (
    fox.pilot$subjectivethreat == 4 ~ 1, 
    fox.pilot$subjectivethreat == 3 ~ .75, 
    fox.pilot$subjectivethreat == 2 ~ .5, 
    fox.pilot$subjectivethreat == 1 ~ .25, 
    fox.pilot$subjectivethreat == 0 ~ 0
  ))

##Appendix C.8
fox.trust.models <- list (
  "Fox Trust" =  lm(fox.trust ~ admissions, data = fox.pilot), 
  "Fox Trust" = lm_lin(fox.trust ~ admissions, ~ white + college + male + age + fox.viewer
                   + media.trust + pid.scale, data = fox.pilot), 
  "ID Threat" = lm(id.threat ~ admissions, data = fox.pilot), 
  "ID Threat" = lm_lin(id.threat ~ admissions, ~ white + college + male + age + fox.viewer
                   + media.trust + pid.scale, data = fox.pilot))

modelsummary(fox.trust.models,
             stars = TRUE,
             coef_rename = c("admissions" = "Admissions", "white" = "white", 
                             "college" = "College educated", "male" = "male", 
                             "age" = "age", "fox.viewer" = "Fox viewer",
                             "media.trust" = "Media trust", "pid.scale" = "PID Scale"),
             notes = list ('Robust Standard Errors in parentheses',
                           'Source: CloudResearch, Nov. 25th-26th, 2023.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

####Main Analysis: Effect of Treatment on future media choice####
##Modeling effect of treatment on future media choice--multinomial logistic regression
##Base/reference category are traditional media--choices are Fox, MSNBC and CNN,  
#Alternative right-wing source (OANN/Newsmax), or no media 

##Collapse media choice into four bins: 
fox.pilot <- fox.pilot %>%
  mutate(future.choice = case_when (
    fox.pilot$futurechoice == 1 ~ 'Fox', 
    fox.pilot$futurechoice == 2 ~ 'Alt Right', 
    fox.pilot$futurechoice == 3 ~ 'Alt Right', 
    fox.pilot$futurechoice == 4 ~ 'Left', 
    fox.pilot$futurechoice == 9 ~ 'Left',
    fox.pilot$futurechoice == 5 ~ 'Traditional', 
    fox.pilot$futurechoice == 8 ~ 'Traditional', 
    fox.pilot$futurechoice == 10 ~ 'None'
  ))

fox.pilot$future.choice <- factor(fox.pilot$future.choice, ordered = FALSE)

##Set "Traditional" as the base category
fox.pilot$future.choice <- relevel(fox.pilot$future.choice, ref = "Traditional")

##Run simple model (no controls)
##Appendix C.9
multinom.media <- multinom(future.choice ~ admissions, data = fox.pilot,
                           Hess = TRUE)

stargazer(multinom.media, title = "Multinomial Logit Model Results", 
          nobs = TRUE, type = "latex")


##Run model with controls##
##Appendix C.10
multinom.media.covar <- multinom(future.choice ~ admissions + white + college + male + age + fox.viewer
                                 + media.trust + pid.scale, data = fox.pilot,
                           Hess = TRUE)

coef_rename <- c("admissions" = "Admissions", 
                 "white" = "White", 
                 "college" = "College", 
                 "male" = "Male", 
                 "age" = "Age", 
                 "fox.viewer" = "Fox Viewer",
                 "media.trust" = "Media Trust",
                 "pid.scale" = "PID Scale")

# Generate the summary table with stargazer
stargazer(multinom.media.covar, title = "Multinomial Logit Model Results", 
          nobs = TRUE, type = "latex")

####Pilot Test: Testing effectiveness/equivalence of each video####
##Appendix C.5
fox_pilot_filtered <- fox.pilot %>%
  filter(admissions == 0)

models_simple_rudy <- list(
  "Big Lie" = lm(big.lie ~ rudy, data = fox_pilot_filtered),
  "Trump Won" = lm(trump.won ~ rudy, data = fox_pilot_filtered),
  "Dominion Changed" = lm(dominion.changed ~ rudy, data = fox_pilot_filtered),
  "Immigrants Voted" = lm(imm.voted ~ rudy, data = fox_pilot_filtered),
  "Venezuela" = lm(venezuela.company ~ rudy, data = fox_pilot_filtered))

modelsummary(models_simple_rudy, vcov = "robust",
             stars = TRUE,
             coef_map = list(rudy = "Giuliani Video", "(Intercept)" = "Intercept"),
             notes = list('Robust Standard Errors in parentheses',
                          'Source: CloudResearch, Nov. 25th-26th, 2023.'),
             gof_omit = 'DF|Deviance|R2|Log.Lik.|BIC|Std. Errors',
             output = "latex")

####Power Analysis: Main Effect####
##Appendix E

#Set up parameters
alpha <- 0.05  # Significance level
power_levels <- seq(0.95, 0.85, by = -0.01)  # Range of power levels
effect_sizes <- seq(0.1, 0.3, by = 0.02)  # Range of effect sizes

#Create a matrix to store sample size results
sample_size_matrix <- matrix(NA, nrow = length(effect_sizes), ncol = length(power_levels))

#Perform power analysis for different effect sizes and power levels
for (i in seq_along(effect_sizes)) {
  for (j in seq_along(power_levels)) {
    sample_size_matrix[i, j] <- pwr.t.test(d = effect_sizes[i], 
                                           sig.level = alpha, power = power_levels[j], alternative = "two.sided")$n
  }
}

#Convert the matrix to a data frame for visualization
result_df <- as.data.frame(cbind(effect_size = rep(effect_sizes, each = length(power_levels)), 
                                 power = rep(power_levels, length(effect_sizes)), 
                                 sample_size = as.vector(sample_size_matrix)))

#Display the result table
print(result_df)

#Plot the results with reversed legend and adjusted Y-axis ticks

plot <- ggplot(result_df, aes(x = effect_size, y = sample_size, group = as.factor(power), color = as.factor(power))) +
  geom_line() +
  geom_point() +
  labs(title = "Power Analysis for Main Effect",
       x = "Effect Size (Cohen's d)",
       y = "Total Sample Size",
       color = "Power") +
  theme_minimal() +
  scale_color_discrete(name = "Power Level", breaks = rev(levels(as.factor(result_df$power)))) +
  scale_y_continuous(breaks = seq(0, max(result_df$sample_size), by = 500))

# Save the plot as a high-quality image
ggsave("power_analysis_plot.png", plot, width = 8, height = 6, units = "in", 
       dpi = 1000)


